# Climate Policy Analysis - Models and Visualizations
# Clean code for analyses after data import

# Load required libraries
library(plm)
library(did)
library(dplyr)
library(ggplot2)
library(broom)
library(stargazer)
library(tidyr)

# Clear environment
rm(list = ls())

#===============================================================================
# SECTION 1: DATA IMPORT AND SETUP
#===============================================================================

# Import the datasets
pdata <- read.csv("pdata.csv")
diddata <- read.csv("diddata.csv")

# Convert back to panel data frame for plm models
pdata <- pdata.frame(pdata, index = c("stcty", "year"))

# Prepare DID data
diddata <- diddata %>%
  mutate(id = as.numeric(factor(stcty)))

cat("Data loaded successfully!\n")
cat("Panel data dimensions:", nrow(pdata), "observations,", length(unique(pdata$stcty)), "counties\n")
cat("DID data dimensions:", nrow(diddata), "observations\n")

#===============================================================================
# SECTION 2: MODEL VARIABLE DEFINITIONS
#===============================================================================

# Define variable groups for systematic analysis
outcome_var <- "RPL_THEMES"  # Lower values = higher vulnerability

# Policy variables by category
policy_vars <- list(
  total = "total_policy",
  
  # Specific features of interest (from your original analysis)
  features = c("FEATURE_economicresilience",
               "FEATURE_ecosystemsnaturalresources", 
               "FEATURE_governmentbylawsordinancescodes",
               "FEATURE_infrastructurebuilt",
               "FEATURE_socialenvironmentaljustice"),
  
  # Plan types of interest
  plans = c("PLAN_TYPE_adaptationplan",
            "PLAN_TYPE_casestudy", 
            "PLAN_TYPE_climatemitigationdocument",
            "PLAN_TYPE_disasterrecoveryplan",
            "PLAN_TYPE_resilienceplan"),
  
  # Community implementation levels
  community = c("COMMUNITY_town", "COMMUNITY_organization", "COMMUNITY_tribe",
                "COMMUNITY_state", "COMMUNITY_watershed"),
  
  # Impact focus areas
  impacts = c("IMPACTS_extremeheat", "IMPACTS_flooding", "IMPACTS_saltwaterintrusion",
              "IMPACTS_sealevelrise", "IMPACTS_stormsurge")
)

# Control variables
control_vars <- c(
  # PCA components
  "PC1_imp", "PC2_imp", "PC3_imp", "PC4_imp",
  "PC1_plan", "PC2_plan", "PC3_plan", "PC4_plan", "PC5_plan",
  "PC1_fea", "PC2_fea", "PC3_fea", "PC4_fea",
  
  # Socioeconomic controls
  "EPL_POV150", "EPL_UNEMP", "M_HBURD_POP", "M_MINRTY_POP",
  "M_SNGPNT_POP", "M_MOBILE_POP", "M_LIMENG_POP",
  
  # Climate exposure
  "crsi_NaturalEnvironment"
)

#===============================================================================
# SECTION 3: MAIN PANEL DATA MODELS
#===============================================================================

# Function to run systematic model analysis
run_policy_models <- function(policy_group, model_type = "random") {
  
  # Create base control formula
  base_controls <- paste(control_vars, collapse = " + ")
  
  results <- list()
  
  for (policy_var in policy_group) {
    # Skip if variable doesn't exist
    if (!policy_var %in% names(pdata)) {
      cat("Warning: Variable", policy_var, "not found in data\n")
      next
    }
    
    # Create formula
    formula_str <- paste(outcome_var, "~", policy_var, "+", base_controls)
    formula_obj <- as.formula(formula_str)
    
    # Run model
    if (model_type == "random") {
      model <- plm(formula_obj, data = pdata, model = "random", random.method = "swar")
    } else {
      model <- plm(formula_obj, data = pdata, model = "within")
    }
    
    # Store results
    results[[policy_var]] <- model
  }
  
  return(results)
}

# Run models for each policy category
cat("\n=== RUNNING PANEL DATA MODELS ===\n")

# 1. Overall policy effect
total_model <- plm(as.formula(paste(outcome_var, "~", policy_vars$total, "+", 
                                    paste(control_vars, collapse = " + "))),
                   data = pdata, model = "random", random.method = "swar")

cat("Total policy model completed\n")

# 2. Policy features
feature_models <- run_policy_models(policy_vars$features)
cat("Feature models completed:", length(feature_models), "models\n")

# 3. Plan types  
plan_models <- run_policy_models(policy_vars$plans)
cat("Plan type models completed:", length(plan_models), "models\n")

# 4. Community levels
community_model <- plm(as.formula(paste(outcome_var, "~", 
                                        paste(policy_vars$community, collapse = " + "), "+",
                                        paste(control_vars, collapse = " + "))),
                       data = pdata, model = "random", random.method = "swar")
cat("Community model completed\n")

# 5. Impact focus areas
impact_models <- run_policy_models(policy_vars$impacts)
cat("Impact models completed:", length(impact_models), "models\n")

#===============================================================================
# SECTION 4: RESULTS EXTRACTION AND VISUALIZATION
#===============================================================================

# Function to extract coefficients for plotting
extract_coefficients <- function(model_list, var_pattern) {
  
  coef_data <- map_dfr(names(model_list), function(var_name) {
    model <- model_list[[var_name]]
    coef_table <- tidy(model)
    
    # Extract coefficient for the policy variable
    policy_coef <- coef_table %>% 
      filter(str_detect(term, var_pattern)) %>%
      slice(1)  # Take first match if multiple
    
    if (nrow(policy_coef) > 0) {
      return(data.frame(
        variable = var_name,
        estimate = policy_coef$estimate,
        std_error = policy_coef$std.error,
        p_value = policy_coef$p.value
      ))
    }
    return(NULL)
  })
  
  return(coef_data)
}

# Extract coefficients for each model category
cat("\n=== EXTRACTING RESULTS ===\n")

# Features
feature_coefs <- extract_coefficients(feature_models, "FEATURE_")
feature_coefs$category <- "Features"

# Plans  
plan_coefs <- extract_coefficients(plan_models, "PLAN_TYPE_")
plan_coefs$category <- "Plan Types"

# Impacts
impact_coefs <- extract_coefficients(impact_models, "IMPACTS_")
impact_coefs$category <- "Impact Focus"

# Community (extract from single model)
community_coefs <- tidy(community_model) %>%
  filter(str_detect(term, "COMMUNITY_")) %>%
  dplyr::select(variable = term, estimate, std_error = std.error, p_value = p.value) %>%
  mutate(category = "Implementation Level")

# Combine all results
all_coefs <- bind_rows(feature_coefs, plan_coefs, impact_coefs, community_coefs)

#===============================================================================
# SECTION 5: VISUALIZATION FUNCTIONS
#===============================================================================

# Function to create coefficient plots
create_coef_plot <- function(coef_data, title, colors = NULL) {
  
  # Clean variable names for display
  coef_data <- coef_data %>%
    mutate(
      clean_name = str_remove(variable, "^(FEATURE_|PLAN_TYPE_|IMPACTS_|COMMUNITY_)"),
      clean_name = str_replace_all(clean_name, c(
        "governmentbylawsordinancescodes" = "Gov. Laws & Codes",
        "infrastructurebuilt" = "Infrastructure Built", 
        "socialenvironmentaljustice" = "Environmental Justice",
        "economicresilience" = "Economic Resilience",
        "ecosystemsnaturalresources" = "Natural Resources",
        "adaptationplan" = "Adaptation Plan",
        "climatemitigationdocument" = "Climate Mitigation",
        "disasterrecoveryplan" = "Disaster Recovery", 
        "resilienceplan" = "Resilience Plan",
        "casestudy" = "Case Study",
        "extremeheat" = "Extreme Heat",
        "saltwaterintrusion" = "Saltwater Intrusion",
        "sealevelrise" = "Sea Level Rise",
        "stormsurge" = "Storm Surge",
        "flooding" = "Flooding",
        "town" = "Town",
        "organization" = "Organization", 
        "tribe" = "Tribe",
        "state" = "State",
        "watershed" = "Watershed"
      )),
      significant = p_value < 0.05
    )
  
  # Default colors if not provided
  if (is.null(colors)) {
    colors <- c("#17a2b8", "#76d7c4", "#20c997", "#ffdd44", "#ffcc00")
  }
  
  # Create plot
  p <- ggplot(coef_data, aes(x = reorder(clean_name, estimate), y = estimate, 
                             color = clean_name, shape = significant)) +
    geom_point(size = 3) +
    geom_errorbar(aes(ymin = estimate - 1.96*std_error, 
                      ymax = estimate + 1.96*std_error), 
                  width = 0.2) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    coord_flip() +
    labs(
      title = title,
      x = "",
      y = "Coefficient Estimate", 
      caption = "Note: Lower RPL_THEMES values indicate higher vulnerability.\nNegative coefficients suggest vulnerability reduction."
    ) +
    scale_shape_manual(values = c(1, 16), name = "Significant", 
                       labels = c("No", "Yes")) +
    theme_minimal(base_family = "Times") +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5),
      legend.position = "bottom",
      axis.text = element_text(size = 10)
    )
  
  # Add colors if multiple categories
  if (length(unique(coef_data$clean_name)) <= length(colors)) {
    p <- p + scale_color_manual(values = colors[1:length(unique(coef_data$clean_name))],
                                guide = "none")
  }
  
  return(p)
}

#===============================================================================
# SECTION 6: GENERATE VISUALIZATIONS  
#===============================================================================

cat("\n=== CREATING VISUALIZATIONS ===\n")

# 1. Policy Features Plot
if (nrow(feature_coefs) > 0) {
  plot_features <- create_coef_plot(
    feature_coefs, 
    "Impact of Policy Features on Climate Vulnerability"
  )
  print(plot_features)
}

# 2. Plan Types Plot
if (nrow(plan_coefs) > 0) {
  plot_plans <- create_coef_plot(
    plan_coefs,
    "Impact of Plan Types on Climate Vulnerability" 
  )
  print(plot_plans)
}

# 3. Impact Focus Plot
if (nrow(impact_coefs) > 0) {
  plot_impacts <- create_coef_plot(
    impact_coefs,
    "Impact of Climate Impact Focus on Vulnerability"
  )
  print(plot_impacts)
}

# 4. Implementation Level Plot
if (nrow(community_coefs) > 0) {
  plot_community <- create_coef_plot(
    community_coefs,
    "Impact of Implementation Level on Climate Vulnerability"
  )
  print(plot_community)
}

# 5. Combined plot
if (nrow(all_coefs) > 0) {
  plot_combined <- ggplot(all_coefs, aes(x = reorder(variable, estimate), 
                                         y = estimate, color = category)) +
    geom_point(size = 2) +
    geom_errorbar(aes(ymin = estimate - 1.96*std_error, 
                      ymax = estimate + 1.96*std_error), 
                  width = 0.2) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    coord_flip() +
    facet_wrap(~category, scales = "free_y") +
    labs(
      title = "Climate Policy Effects on Vulnerability by Category",
      x = "", 
      y = "Coefficient Estimate"
    ) +
    theme_minimal(base_family = "Times") +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5),
      strip.text = element_text(face = "bold"),
      axis.text.y = element_text(size = 8),
      legend.position = "none"
    )
  
  print(plot_combined)
}

#===============================================================================
# SECTION 7: SUMMARY TABLES
#===============================================================================

cat("\n=== GENERATING SUMMARY TABLES ===\n")

# Create summary table of key results
summary_table <- all_coefs %>%
  mutate(
    effect_size = case_when(
      abs(estimate) < 0.01 ~ "Small",
      abs(estimate) < 0.05 ~ "Medium", 
      TRUE ~ "Large"
    ),
    direction = ifelse(estimate < 0, "Reduces Vulnerability", "Increases Vulnerability"),
    significant = ifelse(p_value < 0.05, "Yes", "No")
  ) %>%
  dplyr::select(Category = category, Variable = variable, Estimate = estimate, 
         `Std Error` = std_error, `P-Value` = p_value, 
         `Effect Size` = effect_size, Direction = direction, Significant = significant) %>%
  arrange(Category, Estimate)

print(summary_table)

# Export results
write.csv(summary_table, "policy_effects_summary.csv", row.names = FALSE)
write.csv(all_coefs, "all_coefficients.csv", row.names = FALSE)

#===============================================================================
# SECTION 8: ROBUSTNESS - FIXED EFFECTS MODELS
#===============================================================================

cat("\n=== RUNNING ROBUSTNESS CHECKS (FIXED EFFECTS) ===\n")

# Run key models with fixed effects for robustness
fe_total_model <- plm(as.formula(paste(outcome_var, "~", policy_vars$total, "+", 
                                       paste(control_vars, collapse = " + "))),
                      data = pdata, model = "within")

# Compare random vs fixed effects for total policy effect
comparison_table <- data.frame(
  Model = c("Random Effects", "Fixed Effects"),
  Estimate = c(coef(total_model)["total_policy"], coef(fe_total_model)["total_policy"]),
  Std_Error = c(summary(total_model)$coefficients["total_policy", "Std. Error"],
                summary(fe_total_model)$coefficients["total_policy", "Std. Error"])
)

print("Model Comparison - Total Policy Effect:")
print(comparison_table)

#===============================================================================
# SECTION 9: DID TREATMENT TIMING PREPARATION
#===============================================================================

cat("\n=== PREPARING DID TREATMENT TIMING ===\n")

# Read original RAINE data to get first policy adoption year by county
if (file.exists("pdata.csv")) {
  
  # Create treatment timing for overall policy adoption
  diddata <- diddata %>%
    group_by(stcty) %>%
    arrange(year) %>%
    mutate(
      # First year with any policy
      first_treat = ifelse(any(total_policy > 0), 
                           min(year[total_policy > 0], na.rm = TRUE), 
                           0)
    ) %>%
    ungroup()
  
  # Create treatment timing for specific policy features
  specific_policies <- c(
    "FEATURE_economicresilience",
    "FEATURE_infrastructurebuilt", 
    "FEATURE_governmentbylawsordinancescodes",
    "PLAN_TYPE_adaptationplan",
    "PLAN_TYPE_resilienceplan"
  )
  
  for (policy in specific_policies) {
    if (policy %in% names(diddata)) {
      treat_var <- paste0("treat_", policy)
      diddata <- diddata %>%
        group_by(stcty) %>%
        mutate(
          !!treat_var := ifelse(any(get(policy) > 0, na.rm = TRUE),
                                min(year[get(policy) > 0], na.rm = TRUE),
                                0)
        ) %>%
        ungroup()
    }
  }
  
  # Ensure ID variable exists
  if (!"id" %in% names(diddata)) {
    diddata <- diddata %>%
      mutate(id = as.numeric(factor(stcty)))
  }
  
  cat("Treatment timing variables created\n")
}

#===============================================================================
# SECTION 10: CALLAWAY-SANT'ANNA DID ANALYSIS
#===============================================================================

cat("\n=== RUNNING CALLAWAY-SANT'ANNA DID ANALYSIS ===\n")

# Function to run DID analysis with error handling
run_did_analysis <- function(data, treatment_var, outcome = "SPL_THEMES", 
                             title_suffix = "") {
  
  tryCatch({
    # Check if treatment variable exists and has variation
    if (!treatment_var %in% names(data)) {
      cat("Warning: Treatment variable", treatment_var, "not found\n")
      return(NULL)
    }
    
    treat_summary <- table(data[[treatment_var]])
    if (length(treat_summary) < 2) {
      cat("Warning: No treatment variation for", treatment_var, "\n")
      return(NULL)
    }
    
    cat("Running DID for:", treatment_var, "\n")
    cat("Treatment distribution:", paste(names(treat_summary), "=", treat_summary, collapse = ", "), "\n")
    
    # Run att_gt
    att_result <- att_gt(
      yname = outcome,
      tname = "year", 
      idname = "id",
      gname = treatment_var,
      data = data,
      control_group = "notyettreated",
      bstrap = TRUE,
      panel = TRUE,
      allow_unbalanced_panel = TRUE
    )
    
    # Simple ATT
    simple_att <- aggte(att_result, type = "simple", na.rm = TRUE)
    
    # Dynamic effects (event study)
    dynamic_att <- aggte(att_result, type = "dynamic", na.rm = TRUE)
    
    return(list(
      att_result = att_result,
      simple = simple_att, 
      dynamic = dynamic_att,
      treatment_var = treatment_var
    ))
    
  }, error = function(e) {
    cat("Error in DID analysis for", treatment_var, ":", e$message, "\n")
    return(NULL)
  })
}

# Function to create DID event study plot
create_did_plot <- function(did_result, max_periods = 10) {
  
  if (is.null(did_result) || is.null(did_result$dynamic)) return(NULL)
  
  # Filter to reasonable time window
  event_times <- did_result$dynamic$egt
  keep_periods <- abs(event_times) <= max_periods
  
  if (sum(keep_periods) == 0) return(NULL)
  
  # Create filtered result object
  filtered_result <- list(
    overall.att = did_result$dynamic$overall.att,
    overall.se = did_result$dynamic$overall.se,
    type = did_result$dynamic$type,
    egt = event_times[keep_periods],
    att.egt = did_result$dynamic$att.egt[keep_periods],
    se.egt = did_result$dynamic$se.egt[keep_periods],
    crit.val.egt = did_result$dynamic$crit.val.egt,
    inf.function = did_result$dynamic$inf.function,
    min_e = did_result$dynamic$min_e,
    max_e = did_result$dynamic$max_e,
    balance_e = did_result$dynamic$balance_e,
    call = did_result$dynamic$call,
    DIDparams = did_result$dynamic$DIDparams
  )
  class(filtered_result) <- "AGGTEobj"
  
  # Create clean variable name for title
  clean_name <- str_remove(did_result$treatment_var, "^(treat_|FEATURE_|PLAN_TYPE_)")
  clean_name <- str_replace_all(clean_name, c(
    "governmentbylawsordinancescodes" = "Government Laws & Codes",
    "infrastructurebuilt" = "Infrastructure Built",
    "economicresilience" = "Economic Resilience", 
    "adaptationplan" = "Adaptation Plans",
    "resilienceplan" = "Resilience Plans"
  ))
  
  # Create plot
  plot <- ggdid(filtered_result) +
    labs(
      title = paste("Event Study:", clean_name),
      x = "Years Relative to Policy Adoption",
      y = "Average Treatment Effect",
      caption = "95% confidence intervals shown"
    ) +
    theme_minimal(base_family = "Times") +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5),
      panel.grid.minor = element_blank()
    ) +
    scale_color_manual(
      values = c("#ffdd44", "#76d7c4"),
      labels = c("Pre-Treatment", "Post-Treatment")
    )
  
  return(plot)
}

# 1. Overall policy adoption DID
cat("\n--- Overall Policy Adoption ---\n")
overall_did <- run_did_analysis(diddata, "first_treat")

if (!is.null(overall_did)) {
  cat("Simple ATT:\n")
  print(summary(overall_did$simple))
  
  # Create and display plot
  overall_plot <- create_did_plot(overall_did)
  if (!is.null(overall_plot)) {
    print(overall_plot)
  }
}

# 2. Specific policy feature DID analyses
did_results <- list()

specific_treatments <- paste0("treat_", specific_policies)
names(specific_treatments) <- specific_policies

for (i in seq_along(specific_treatments)) {
  policy_name <- names(specific_treatments)[i]
  treat_var <- specific_treatments[i]
  
  cat("\n---", policy_name, "---\n")
  
  did_result <- run_did_analysis(diddata, treat_var)
  
  if (!is.null(did_result)) {
    did_results[[policy_name]] <- did_result
    
    # Print simple ATT
    cat("Simple ATT:\n")
    print(summary(did_result$simple))
    
    # Create and display plot
    did_plot <- create_did_plot(did_result)
    if (!is.null(did_plot)) {
      print(did_plot)
    }
  }
}

#===============================================================================
# SECTION 11: DID RESULTS SUMMARY
#===============================================================================

cat("\n=== DID RESULTS SUMMARY ===\n")

# Extract simple ATT results
extract_did_summary <- function(did_result, name) {
  if (is.null(did_result) || is.null(did_result$simple)) return(NULL)
  
  summary_obj <- summary(did_result$simple)
  
  data.frame(
    policy = name,
    att = summary_obj$overall.att,
    se = summary_obj$overall.se, 
    pvalue = 2 * (1 - pnorm(abs(summary_obj$overall.att / summary_obj$overall.se))),
    n_treated = length(unique(did_result$att_result$group[did_result$att_result$group > 0])),
    stringsAsFactors = FALSE
  )
}

# Compile DID summary table
did_summary_list <- list()

# Overall results
if (!is.null(overall_did)) {
  did_summary_list[["Overall Policy"]] <- extract_did_summary(overall_did, "Overall Policy")
}

# Specific policy results  
for (policy_name in names(did_results)) {
  if (!is.null(did_results[[policy_name]])) {
    clean_name <- str_remove(policy_name, "^(FEATURE_|PLAN_TYPE_)")
    clean_name <- str_to_title(str_replace_all(clean_name, c(
      "governmentbylawsordinancescodes" = "gov_laws_codes",
      "infrastructurebuilt" = "infrastructure", 
      "economicresilience" = "economic_resilience",
      "adaptationplan" = "adaptation_plan",
      "resilienceplan" = "resilience_plan"
    )))
    
    did_summary_list[[clean_name]] <- extract_did_summary(did_results[[policy_name]], clean_name)
  }
}

# Combine results
if (length(did_summary_list) > 0) {
  did_summary_table <- bind_rows(did_summary_list) %>%
    mutate(
      ci_lower = att - 1.96 * se,
      ci_upper = att + 1.96 * se,
      significant = pvalue < 0.05,
      effect_direction = ifelse(att < 0, "Reduces Vulnerability", "Increases Vulnerability")
    ) %>%
    arrange(att)
  
  print("DID Analysis Summary:")
  print(did_summary_table %>% 
          dplyr::select(Policy = policy, ATT = att, `Std Error` = se, 
                 `P-Value` = pvalue, `95% CI Lower` = ci_lower, 
                 `95% CI Upper` = ci_upper, `Effect Direction` = effect_direction))
  
  # Export DID results
  write.csv(did_summary_table, "did_results_summary.csv", row.names = FALSE)
  cat("\nDID results exported to: did_results_summary.csv\n")
}

#===============================================================================
# SECTION 12: MODEL COMPARISON VISUALIZATION
#===============================================================================

cat("\n=== CREATING MODEL COMPARISON PLOTS ===\n")

# Compare panel data vs DID results for overlapping policies
if (exists("did_summary_table") && nrow(all_coefs) > 0) {
  
  # Prepare comparison data
  panel_comparison <- all_coefs %>%
    filter(str_detect(variable, paste(specific_policies, collapse = "|"))) %>%
    mutate(
      clean_name = str_remove(variable, "^(FEATURE_|PLAN_TYPE_)"),
      clean_name = str_to_title(str_replace_all(clean_name, c(
        "governmentbylawsordinancescodes" = "gov_laws_codes",
        "infrastructurebuilt" = "infrastructure",
        "economicresilience" = "economic_resilience", 
        "adaptationplan" = "adaptation_plan",
        "resilienceplan" = "resilience_plan"
      ))),
      method = "Panel Data",
      coefficient = estimate,
      std_error_coef = std_error
    ) %>%
    dplyr::select(policy = clean_name, method, coefficient, std_error_coef)
  
  did_comparison <- did_summary_table %>%
    mutate(
      method = "DID (Callaway-Sant'Anna)",
      coefficient = att,
      std_error_coef = se
    ) %>%
    dplyr::select(policy, method, coefficient, std_error_coef)
  
  # Combine for comparison plot
  comparison_data <- bind_rows(panel_comparison, did_comparison) %>%
    filter(!is.na(coefficient))
  
  if (nrow(comparison_data) > 0) {
    comparison_plot <- ggplot(comparison_data, aes(x = policy, y = coefficient, 
                                                   color = method, shape = method)) +
      geom_point(size = 3, position = position_dodge(width = 0.3)) +
      geom_errorbar(aes(ymin = coefficient - 1.96 * std_error_coef,
                        ymax = coefficient + 1.96 * std_error_coef),
                    width = 0.2, position = position_dodge(width = 0.3)) +
      geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
      coord_flip() +
      labs(
        title = "Model Comparison: Panel Data vs. Difference-in-Differences",
        x = "Policy Type",
        y = "Coefficient Estimate", 
        color = "Method",
        shape = "Method",
        caption = "95% confidence intervals shown"
      ) +
      scale_color_manual(values = c("#e74c3c", "#3498db")) +
      theme_minimal(base_family = "Times") +
      theme(
        plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "bottom"
      )
    
    print(comparison_plot)
    
    # Export comparison data
    write.csv(comparison_data, "method_comparison.csv", row.names = FALSE)
  }
}

